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We consider a class of neural circuit models with internal noise sources arising in sensory 
systems. The basic neuron model in these circuits consists of a dendritic stimulus 
processor (DSP) cascaded with a biophysical spike generator (BSG). The dendritic stimulus 
processor is modeled as a set of nonlinear operators that are assumed to have a Volterra 
series representation. Biophysical point neuron models, such as the Hodgkin-Huxley 
neuron, are used to model the spike generator. We address the question of how intrinsic 
noise sources affect the precision in encoding and decoding of sensory stimuli and the 
functional identification of its sensory circuits. We investigate two intrinsic noise sources 
arising (i) in the active dendritic trees underlying the DSPs, and (ii) in the ion channels 
of the BSGs. Noise in dendritic stimulus processing arises from a combined effect of 
variability in synaptic transmission and dendritic interactions. Channel noise arises in the 
BSGs due to the fluctuation of the number of the active ion channels. Using a stochastic 
differential equations formalism we show that encoding with a neuron model consisting 
of a nonlinear DSP cascaded with a BSG with intrinsic noise sources can be treated 
as generalized sampling with noisy measurements. For single-input multi-output neural 
circuit models with feedforward, feedback and cross-feedback DSPs cascaded with BSGs 
we theoretically analyze the effect of noise sources on stimulus decoding. Building on 
a key duality property, the effect of noise parameters on the precision of the functional 
identification of the complete neural circuit with DSP/BSG neuron models is given. We 
demonstrate through extensive simulations the effects of noise on encoding stimuli with 
circuits that include neuron models that are akin to those commonly seen in sensory 
systems, e.g., complex cells in V1. 



Keywords: Volterra dendritic stimulus processors, biophysical spike generators, noise, neural encoding, neural 
decoding, functional identification, Hodgkin-Huxley neuron, phase response curve 



1. INTRODUCTION 

Intrinsic noise sources are diverse and appear on many levels 
of a neuronal system ranging from electrical to chemical noise 
sources (Faisal et al., 2008; Destexhe and Rudolph-Lilith, 2012) 
and from single cells to networks of neurons. At the cellular 
and subcellular level, variability in biochemical reactions leads 
to stochastic transduction processes (Song et al, 2012), and ion 
channel fluctuations (Neher and Sakmann, 1976; White et al., 
1998) result in variability in spike generation and propagation 
(Faisal and Laughlin, 2007). At the network level, probabilistic 
quantal release of neurotransmitters (Katz, 1962), background 
synaptic activity (Destexhe et al, 2003; Jocobson et al., 2005) and 
variability in timing of spikes from presynaptic neurons (Faisal 
and Neishabouri, 2014) are sources of stochastic fluctuation of 
synaptic conductances (Destexhe et al, 2001) that are believed 
to have a major impact on spike time variability (Yarom and 
Hounsgaard, 2011). 

The existence of sources of noise also leads to variability in the 
spike times even when neurons are subject to the same, repeated 
inputs (Calvin and Stevens, 1968; Berry et al, 1997; de Ruyter van 
Steveninck et al, 1997). Spikes are the primary form of carriers of 



information in the nervous system and their timing is thought to 
be relevant to the message neurons need to convey (Rieke et al., 
1999). Therefore, the variability of spike timing may reduce or 
damage the information being transmitted. It is quite remarkable, 
however, that sensory systems manage to be very robust even if 
they are subject to interference due to noise. Visual and auditory 
systems are two examples in which the stimuli are highly time 
varying. These systems have been reported to convey information 
with high spike timing precision (Butts et al., 2007; Kayser et al., 
2010). 

Noise maybe useful in facilitating signal detection (McDonnell 
and Ward, 2011). Still, interference due to noise poses an impor- 
tant limit on how well sensory systems can represent input 
stimuli. It is not clear how intrinsic noise sources affect the rep- 
resentation of sensory inputs based on spike times, and how they 
impact the functional identification of sensory neurons. 

We study the representation of sensory stimuli using a novel 
neural circuit model, that extends previously proposed mod- 
els (Lazar et al, 2010; Lazar and Slutskiy, 2014, in press) in 
terms of architectural complexity and the existence of intrin- 
sic noise sources. Our base level circuit architecture consists of 
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two interconnected neurons, each with two cascaded stages. The 
first stage comprises two types of dendritic stimulus processors. 
The first dendritic stimulus processor performs nonlinear pro- 
cessing of input stimuli in the feedforward path leading to the 
spike generator. The second dendritic stimulus processor per- 
forms nonlinear processing in the feedback loop whose inputs are 
spike trains generated by biophysical spike generators (BSGs). The 
BSGs constitute the second stage of the base level circuit. 

Our nonlinear dendritic stimulus processors describe func- 
tional I/O relationships between the dendritic outputs in the first 
stage and inputs that are either sensory stimuli or spikes gener- 
ated by BSGs. DSPs are modeled using Volterra series. Volterra 
series have been used for analyzing nonlinear neuronal responses 
in many contexts (Lu et al., 2011; Eikenberry and Marmarelis, 
2012), and have been applied to the identification of single neu- 
rons in many of sensory areas (Benardete and Kaplan, 1997; 
Theunissen et al., 2000; Clark et al., 2011). Volterra dendritic 
processors can model a wide range of nonlinear effects com- 
monly seen in sensory systems (Lazar and Slutskiy, in press). 
Here, in addition, we introduce nonlinear interactions between 
neurons in the feedback and cross-feedback paths. This gives 
rise to interesting neural processing capabilities directly in the 
spike domain, e.g., coincidence detection (Agmon-Snir et al., 
1998; Stuart and Hausser, 2001). The relationships described 
here by the Volterra model are functional and do not address 
the underlying circuit/dendritic tree level interactions. However, 
the latter have recently been subject to intense investigations 
(London and Hausser, 2005; Wohrer and Kornprobst, 2009; 
Werblin, 2011; Xu et al, 2012; Yonehara et al, 2013; Zhang 
et al, 2013). Conductance-based, biophysical spike generators are 
well established models that have been extensively used in stud- 
ies of neuronal excitability and in large simulations of spiking 
neural networks (Izhikevich, 2007). Following Lazar (2010), we 
use formal BSG models to represent sensory stimuli under noisy 
conditions. 

We formulate the encoding, decoding and functional iden- 
tification problems under the neural encoding framework of 
Time Encoding Machines (TEMs). In this modeling framework 
the exact timing of spikes is considered to carry information 
about input stimuli (Lazar and Toth, 2004). The separation into 
dendritic stimulus processors and spike mechanisms mentioned 
above allows us to study synaptic inputs and spike genera- 
tion mechanisms separately, and hence independently model the 
intrinsic noise sources of each component. We incorporate two 
important noise sources into a general single-input multi-output 
neural circuit model. The first is a channel noise source that arises 
in spike generation (White et al, 2000). The second is a synap- 
tic noise source due to a variety of fluctuating synaptic currents 
(Manwani and Koch, 1999). 

Based on the rigorous formalism of TEMs, we show how noise 
arising in dendritic stimulus processors and in biophysical spike 
generators is related to the measurement error in generalized sam- 
pling. Dendritic stimulus processing and spike generation can 
then be viewed as a generalized sampling scheme that neurons 
utilize to represent sensory inputs (Lazar et al., 2010). Contrary 
to traditional sampling where the signal amplitude is sampled at 
clock times, neurons asynchronously sample all stimuli. 



We systematically investigate how the strength of noise sources 
degrades the faithfulness of stimulus representation and the 
quality of functional identification of our proposed class of neural 
circuits. Furthermore, since the representation is based on spike 
timing, it is natural to investigate how spike timing variability 
affects the precision in representing the amplitude information 
of sensory stimuli. 

The work presented here requires a substantial amount of 
investment in the mathematical formalism employed throughout. 
There are a number of benefits in doing so, however. Formulating 
the problem of stimulus encoding with a neural circuit with 
intrinsic noise sources as one of generalized sampling, i.e., of tak- 
ing noisy measurements is of interest to both experimentalists 
and theoreticians alike. Understanding that the problem of neu- 
ral decoding and functional identification are dual to each other 
is key to building on either or both. Finding how many repeat 
experiments need to be performed for a precise quantitative iden- 
tification of Volterra kernels is of great value in neurophysiology. 
A further qualitative insight of our work is that for neural cir- 
cuits with arbitrary connectivity, feedforward kernels are typically 
easier to estimate than feedback kernels. Finally, our finding that 
some key nonlinear neural circuits are tractable for detailed noise 
analysis suggests a wide reaching analytical methodology. 

2. MODELING NONLINEAR NEURAL CIRCUITS, STIMULI, 
AND NOISE 

We present in Section 2.1 the general architecture of the neural 
circuits considered in this paper. In Section 2.2 we discuss the 
modeling of the space of stimuli. Volterra DSPs are the object 
of Section 2.3. Finally, in Section 2.4 we provide models of BSGs 
with intrinsic noise sources. 

2.1. NEURAL CIRCUIT ARCHITECTURE 

The general architecture of the neural circuit considered here is 
shown in simplified form in Figure 1 . It consists of two neu- 
rons with a common time-varying input stimulus. With added 
notational complexity the neural circuit in Figure 1 can easily be 
extended in two ways. First, multiples of such circuits can encode 
a stimulus in parallel (see Section 2.1 in the Supplementary 
Material). In this case only pairs of neurons are interconnected 
through the feedback kernels. Second, more neurons can be con- 
sidered in the neural circuit of Figure 1; all these neurons can be 
fully interconnected through feedback loops. 

Each neuron i, i = 1,2, receives a single time-varying input 
stimulus «i (f). The modeling of the input stimulus is discussed in 
Section 2.2. The output of each of the biophysical spike generators 
(BSGs) is a spike sequence denoted by (tj) and (f ( 2 ), k, I € Z. 

The input stimulus ui(t) is first processed by a feedfor- 
ward Dendritic Stimulus Processor (feedforward DSP) (Lazar 
and Slutskiy, in press). The feedforward DSP models the aggre- 
gated effect of processing in the neural circuits in the prior 
stages and in the dendritic tree of neuron i = 1,2. For exam- 
ple, if the neurons in the model circuit are considered to be 
Retinal Ganglion Cells (RGCs), then the feedforward Volterra 
DSP models the processing that takes place in the outer- and 
inner-plexiform layers of the retina as well as in the dendritic 
trees of an RGC (Werblin, 201 1; Masland, 2012). The feedforward 
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FIGURE 1 | Diagram of the architecture of the neural circuits. 



DSPs are modeled here as second order Volterra expansion terms 
(Volterra, 1930). The first order terms h\ u (t) in the feedforward 
DSPs are linear filters typically used in modeling receptive fields. 
The second order terms hi 1 * (fi, h) model nonlinear operations 
on the stimulus uy(t). 

A second group of Volterra DSPs models the cross-feedback 
interactions between the two neurons. Instead of time-varying 
stimuli, the output spikes generated by the BSGs are the inputs 
to these DSPs. We therefore refer to these as feedback Dendritic 
Stimulus Processors (feedback DSPs). The output spikes of 
each individual neuron i are processed by the first order term 
]if(t), i,j= 1,2, J ^= j. In addition, output spikes from both 
neurons interact nonlinearly through the second order terms 
hj'iti, tz), i,j = 1, 2, i ^ j. The summed responses from the 
first order feedback DSP and the second order feedback DSP 

hj' are fed back to neuron i as additional dendritic currents. 

The dendritic currents consisting of the output of the DSPs 
with added noise are subsequently encoded by biophysical spike 
generators. BSGs are biophysically realistic axon hillock spike gen- 
erator models that are governed by a set of differential equations 
with multiple types of ion channels (Hodgkin and Huxley, 1952; 
Izhikevich, 2007). The detailed BSG models are introduced in 
Section 2.4. The spike times of output spikes generated by the 
BSGs are assumed to be observable. 

We identify two intrinsic noise sources of the proposed neu- 
ral circuit. First, the feedforward DSPs and the feedback DSPs are 
affected by additive Gaussian white noise. This noise arises from 
the combined effect along the path from sensory transduction 
to synaptic integration and includes synaptic background noise 
and stochasticity in the dendritic tree (Manwani and Koch, 1999; 



Fellous et al., 2003; Destexhe and Rudolph-Lilith, 2012). Since 
the outputs of the feedforward and feedback DSPs are additively 
combined, we consider, for simplicity, a single source of additive 
Gaussian white noise. Second, the ion channels of the BSGs are 
intrinsically stochastic and introduce noise in the spike generators 
(White et al, 2000; Hille, 2001). 

2.2. MODELING SIGNAL SPACES 

Two signal spaces will be considered here. The first, models the 
space of input signals to feedforward DSPs. The second models 
the space of input spikes to feedback DSPs. These spaces will be 
formally described below. 

2.2. 1. Modeling the space of input stimuli 

We model the space of input stimuli as a Reproducing Kernel 
Hilbert Space (RKHS) (Berlinet and Thomas-Agnan, 2004). 
RKHSs are versatile vector spaces for modeling signals arising 
in computational neuroscience, signal processing and machine 
learning. For example, auditory signals, olfactory signals and 
visual signals can readily be modeled as band-limited functions of 
an RKHS with a sine or Dirichlet kernel (Lazar et al., 2010; Lazar 
and Slutskiy, 2013). A particular choice of RKHSs in this arti- 
cle is the space of trigonometric polynomials. The computational 
advantage of working on the space of trignometric polynomi- 
als has been discussed (Lazar et al., 2010) and is closely related 
to the algorithmic tractability of the Fourier series in the dig- 
ital domain. If the biological signals have unknown bandwidth 
with a spectrum that falls off fast enough, many Sobolev spaces 
might be a suitable choice of RKHS (Berlinet and Thomas- 
Agnan, 2004; Lazar and Pnevmatikakis, 2009). In such spaces the 
norm may include the derivative of the signal, i.e., the rate of 
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change of the signal that many neurons are sensitive to Kim et al. 
(2011). 

The space of trigonometric polynomials is defined as below. 

Definition 2.1. The space of trigonometric polynomials 7i\ is a 
function space whose elements are functions defined on the domain 
Di = [0, S 1 ], S l e R+, of the form 

i 1 

ui(t) = uiei{t), (1) 
X=-i> 

•where 

er(t) = -^=ffi\l= -I 1 ,-.- .I 1 , (2) 

are a set of orthonormal basis functions. S2 1 denotes the bandwidth 
and L 1 is the order of the space. 

7i\ endowed with the inner product: 



(Uy,vi) = / ui(t)vi{t)dt (3) 

is a Hilbert Space. Intuitively, the basis functions £/( f )> 
I = —L 1 , ... ,L , can be interpreted as a set of discrete spectral 
lines uniformly spaced in the frequency domain between — Q} 
and Q, 1 . For a given signal U\ (t), the amplitude of its spectral lines 
is determined by the coefficients U\,l= —L 1 , ,..,!}. 

Remark 2.2. Functions in 7i\ are periodic over K with period 

S 1 = ^f~- Therefore, the domain Di covers exactly one period of 
the function. Note that the ufs are closely related to the Fourier 
coefficients of the periodic signal u\(t), and can thereby be very 
efficiently computed via the Fast Fourier Transform. 

7i\ is an RKHS with reproducing kernel (RK) 

I 1 

;=-!' 

It can be easily verified that the RK satisfies the reproducing 
property 

(ui( •),•*! ft •)> = «iW,Vui eH\,te Di. (5) 



Definition 2.3. We shall also consider the tensor product space Ti\ 
on the domain D 2 = [0, S 1 ] x [0, S 1 ], whose elements are of the 
form 

I 1 L [ 

U 2 (ti,t 2 )= Yl Yl u hh e hh^l^2), (6) 
/j = -L 1 h = —I 1 



where 

e h , 2 (tut 2 ) = ^e' h ^ tl e jh ^ t \ (7) 
are a set of functions forming an orthonormal basis. 
Ti\ is again an RKHS with RK 

L 1 L 1 

Kl{h,t 2 ;si,s 2 )= e hh(h -s\,t 2 -s 2 ). (8) 

h = -V-h = -l> 

Note that we use the subscript to indicate the dimension of the 
domain of functions, i.e., the number of variables the functions 
in the RKHS have, and use the superscript 1 to indicate the input 
space. 

Projections of functions onto the RKHSs introduced here can 
be defined as follows: 

Definition 2.4. Let h\ e L'(Di), where L 1 denotes the space of 
Lebesgue integrable functions. The operator V 1 : L'(Di) — > Ti.\ 
given by 

(V 1 h 1 )(t)= f h(s)K{(t; s)ds, 

is called the projection operator from L (Ui) to Ji\. Similarly, let 
h 2 (t u t 2 ) e V(B 2 ), the operator V 1 : L 1 (D 2 ) ->• U\ (by abuse 
of notation ) given by 

(V 1 h 2 )te,t 2 )= / h 2 (s l ,s 2 )K\(t 1 ,t 2 ;si,s 2 )ds l ds 2 , 

is called the projection operator from L 1 (D 2 ) to H\. 
2.2.2. Modeling the space of spikes 

The feedback kernels of the neural circuit in Figure 1 receive as 
inputs spike trains generated by the BSGs. Spike trains are often 
modeled as sequences of Dirac delta pulses and, consequently, the 
outputs of linear feedback kernels are the result of superposition 
of their impulse responses (Keat et al., 2001; Pillow et al., 2008; 
Lazar et al, 2010). 

Dirac delta pulses have infinite bandwidth. Spikes generated by 
the BSGs, however, have limited effective bandwidth. Following 
(Lazar and Slutskiy, 2014) spikes are modeled to be the RK of 
an one-dimensional Hilbert space 7Yf at spike time occurrence. 
Here Tij is a space of trigonometric polynomials whose order I 2 , 
period S 2 and bandwidth £2 2 may differ from the input stimulus 
space H\, where Q. 2 shall be larger than the bandwidth assumed 
for the feedback kernel, and S 2 is much larger than the support of 
the feedback kernel (Lazar and Slutskiy, 2014). A spike at time t' k 
of neuron i can then be expressed in functional form as K 2 (t' k ; t), 
where the superscript indicates that the RK belongs to the spike 
input space. 
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Due to the reproducing property, single or pairs of input spikes 
have the property 

h(t) * K 2 (4; t) = £ h(t- s)K 2 (4; s) ds = (p 2 h) (t - t[) 
and 

jf h 2 (t- s u t - s 2 )Kj (ti <j; si, n) dsids 2 = (V 2 h 2 ) (t -t\,t- {j 

for i, j = 1,2, i / j. The operator V 2 is similarly defined to 
V 1 above; it denotes, however, the projection onto the space 
of spikes. Thus, not surprisingly, incoming spikes directly read- 
out the projection of the feedback kernels. By letting L 2 — > 00, 
(V 2 h\^ (f — fjt) shall converge to h it — h) in L 2 norm as the 
RK converges to the sine function and the RKHS becomes the 
space of band-limited signals (Lazar et al., 2010). A more detailed 
analysis is available in Lazar and Slutskiy (2014). This formalism 
will be employed for solving the functional identification problem 
formulated in Section 4.1. 

2.3. VOLTERRA DENDRITIC STIMULUS PROCESSORS 

As mentioned in Section 2.1, two forms of dendritic stimulus 
processing appear in our model. 

2.3. 1. Feedforward Volterra dendritic stimulus processors 

The feedforward DSPs are modeled as up to second order terms in 
the Volterra series. The feedforward DSPs take continuous signals 
in the stimulus space as inputs, while the output can be expressed 
as (see also Figure 1) 



it — s)u\ is)ds + 



/> 

J0 2 



it - si, t - s 2 ) uiisi)uiis 2 )ds\ds 2 , (9) 



where h\ u e L^Di) and h\ u e ViB 2 ) denote, respectively, 
the first and second order Volterra kernels, i = 1,2. They are 
assumed to be real, causal and bounded-input bounded-output 
(BIBO)-stable. It is also assumed that both h\ and have finite 
memory. In addition, h\ l1 is assumed, without loss of generality, 
to be symmetric, i.e., hi (h, t 2 ) = h^ 1 ' (f 2 , h). 

Example 2.5. We present here a Volterra DSP that is akin to a 
model of dendritic stimulus processing of complex cells in the pri- 
mary visual cortex (VI ). The difference is that the complex cells 
operate spatio-temporally, whereas in the example given below they 
operate temporally. We first consider two first order kernels based on 
Gabor functions, 

g c it) = exp (- ^ ~ n °^^ c ) cos(2jr • 10 • (f - 0.13)) , 



gsit) = exp 



2 • 0.0005 

(f-0.13) 2 
2 • 0.0005 



sin(27r • 10 • (f-0.13)) 



The two filters are Gaussian modulated sinusoids, that are typically 
used to model receptive fields of simple cells in the primary visual 
cortex (VI ) where the variables denote space instead of time (Lee, 



1996; Dayan and Abbott, 2001). In addition, the two filters are 
quadrature pair in phase. Both filters are illustrated in Figure 2A. 
The response of applying the input stimulus u\ on the temporal fil- 
ters with impulse response g c and g s is given by f 0i g c (t — s)u\ (s)ds 
and J Di g s (t — s)u\is)ds, respectively. 

The responses of the two linear filters of the complex cell model 
are squared and summed to produce the phase invariant measure v' 
(Carandini et al, 2005), where 



v'it) 



1- 

Jo, 



(t — s)u\(s)ds 



+ 



7 

J'O, 



s (t — s)u\(s)ds 



c (t - si )h (t - s 2 )«i (si )«i (s 2 )dsids 2 



/ 

/ gs 

/ [gcit ~ Si)g c (t - S 2 ) +gs(t - Si)g s (t - S 2 )] 

Jn 2 



(t - si)g s (t - s 2 )ui(si)ui(s 2 )dsids 2 



(10) 



/ 

Js-. 



ui(si)ui(s 2 )dsids 2 



hl u it — si,t— s 2 )ui(si)ui(s 2 )dsids 2 , 



wherehl u {tu t 2 ) = g c (ti)g c {t 2 ) + gs(ti)g s {t 2 ). Therefore, theoper- 
ation performed by a complex cell can be modeled with a second 
order Volterra kernel. h\ h is shown in Figure 2B. 

We now take a closer look at the operation of the second order 
kernel. The two dimensional convolution of the second order kernel 
with u 2 (t\, t 2 ) is shown in Figure 2C. 

It is important to note that, since the second order kernel has 
finite memory, it may not have enough support to cover the entire 
domain H) 2 for u 2 (ti, t 2 ). For example, as illustrated in Figure 2C, 
the output of the second order feedforward DSP at time t is given 
by the integral of the product of u 2 it\, t 2 ) and a rotated h\ l1 with 
the origin shifted to it, t) [see also (10)]. Since the shift is along the 
diagonal, only U 2 (ti, t 2 ) in the domain that is contained within the 
black lines is multipled by nonzero values ofh\ h . u 2 (t\ , t 2 ) elsewhere 
in the domain is always multiplied by zero in evaluating the output. 
Therefore, the output of the second order filter only contains infor- 
mation about u 2 within the domain located in between the black 
lines in Figure 2C. This has implications on decoding the signal (see 
also Remark 3.11 in Section 3.2) 

2.3.2. Feedback Volterra dendritic stimulus processors 

As already mentioned, the feedback DSPs do not operate on stim- 
uli directly but rather on spikes generated by BSGs. We assume 
that/if' € V(Bi), hf € L 1 (D 2 ),i / j, are real, causal, BIBO- 
stable and have finite memory. In addition, we assume that these 
kernels are effectively band-limited (see also Section 2.2.2). In 
functional form we denote a train of spikes as Ylk^-i^k' ^ ne 
output of the feedback DSP i amounts to 



£ (P 2 h 2ji 



)H) 



+ 



E E ^ 



')( 



t - t / , t - 



(ID 



ithj # 
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FIGURE 2 | Examples of Volterra kernels. (A) First order kernels of 
quadrature pair of Gabor functions modeling the receptive fields of simple 
cells. (B) Second order kernel modeling receptive fields of complex cells. (C) 
The mechanics of the two dimensional convolution operation between the u 2 
(S 1 = 0.8, D 2 = [0,0.8] x [0,0.8]) and h^'. u 2 (U,t 2 ) = m(fi)ui(f 2 ) is shown 
in the background. The inset shows the second order Volterra kernel h™' 



rotated 180° around origin [see also (B)]. (hV' is only shown in a restricted 
domain and is zero elsewhere). For f = 0.3, the output of the convolution is 
the integral of the product of the rotated Volterra kernel and the signal 
underneath. Since the convolution is evaluated on the diagonal f = t\ = t 2 , 
the second order kernel shifts, as t increases, along the arrow on the 
diagonal. See also Supplementary Figure 5E. 



In particular, the inputs to the second order term of the 
feedback DSPs are generated by two neurons. This allows for 
modeling nonlinear interactions between the two neurons in the 
spike domain. 

2.3.3. Overall output from DSPs 

The overall inputs (without noise) to the two BSGs in Figure 1 are 
v\t)=\ h\ u (t-s)u l (s)ds+ [ h^ i {t-si,t-s 2 )ui{si) 

Ul (s 2 )d Sl ds 2 + £ (v 2 hf ) (t-$ + £E (r 2 hf) 
(t - tj, f - 4) > and »> y = i, 2, f ^ j. 

(12) 

The system of Equations (12) above functionally describe the 
post-synaptic aggregate currents that are injected into the 
BSG i. 

There are a variety of noise sources to be considered. Synaptic 
variability of feedforward DSPs adds noise sources to the cur- 
rent input to the BSGs. These include thermal noise, synap- 
tic background noise, etc. (Jonston, 1927; Calvin and Stevens, 
1968; Manwani and Koch, 1999; Fellous et al, 2003; Destexhe 
and Rudolph-Lilith, 2012). Feedback DSP kernels may them- 
selves be subject to intrinsic noise sources that may lead to 
variability in the spike generation process. Intrinsic variabil- 
ity of BSG spike times can, e.g., contribute to the variability 
of the aggregate current driving the axon hillock in feedback 
loops. 

Overall, the combined effect of DSP noise sources is mod- 
eled as Gaussian white noise processes that are added to the 
feedforward and feedback DSP outputs. The sum total of sig- 
nal and noise represents the aggregate current input to the 
BSGs (see Figure 1). Formal DSP noise models will be incor- 
porated directly into the BSG model presented in the next 
section. 



2.4. BIOPHYSICAL SPIKE GENERATORS 
2.4.1. BSGs and phase response curves 

We consider biophysically realistic spike generators such as 
the Hodgkin-Huxley, Morris-Lecar, Connor- Stevens neurons 
(Hodgkin and Huxley, 1952; Connor and Stevens, 1971; Morris 
and Lecar, 1981). The class of BSGs can be expressed in vector 
notation as 

f = f( X ^),i=l,2, (13) 

where x 1 are the state variables, f are vector functions of the same 
dimension, and I' are the constant bias currents in the voltage 
equation of each BSG. 

Each input current v'(t) is applied to the neuron i by additive 
coupling to the voltage equation, typically the first of the set of 
ordinary differential equations, i.e., 

^ =f , (x',7') + [v'(f),o] T ,i= 1,2, (14) 

where 0 is a row vector of appropriate size. 

We assume that the neuron is periodically spiking when no 
external input is applied. This can be satisfied by a constant bias 
current I' additively coupled onto the voltage equation. The use 
of P is necessary to formulate the encoding for the single neuron 
case, and this assumption will be relaxed later in this article. 

A large enough bias current induces a periodic oscillation of 
the biophysical spike generator. Therefore, the phase response 
curve (PRC) is well defined for this limit cycle (Izhikevich, 2007). 
We denote the PRC of the limit cycle induced by the bias current 

r as f (t, r) = Ur\ (f, r) , f[ (t, ^ (t, rj\ T with 

appropriate dimension N\ where i/r* (f , 7') , n = 1, 2, • • • , N', 
are the PRCs associated with the nth state variable. Without loss of 
generality, we assume that yjr\ (f, 7 ! ) is always the PRC associated 
with the voltage variable. 

An example of a Hodgkin-Huxley neuron model of a BSG can 
be found in Section 2.2 in the Supplementary Material. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



September 2014 | Volume 8 | Article 95 | 6 



Lazar and Zhou 



Neural circuits with intrinsic noise sources 



, dW^ j , and W p (t),p = 2, ■ ■ ■ , P>, 



2.4.2. Channel noise in BSGs 

As shown in Figure 1, we consider BSGs with noise sources in the 
ion channels. The noise arises due to thermal fluctuations (White 
et al., 2000; Hille, 2001) as the finite number of ion channels in 
the BSGs open and close stochastically. 

The differential equations that govern the dynamics of the 
BSGs in (14) are deterministic. The set of stochastic differen- 
tial equations (SDEs) below represent their stochastic counterpart 
(Lazar, 2010): 

dV = V (Y\ A dt + B' (Y') dZ'(t), i = 1, 2, (15) 

where B' is a matrix with state dependent values, dZ 1 = 
r . . -iT 

Iv'dt, dW'^dW^ 

independent Brownian motion processes. Note that P' does not 
necessarily have to be equal to N', the number of state variables. 
The first element in the stochastic differential dZ' is the aggre- 
gate dendritic input v'dt driving the voltage equation. The other 
entries in dZ' are noise terms that reflect the stochastic fluctuation 
in the ion channels / gating variables. 

Randomness is often added to BSGs by setting B' = I, where 
I is a bP x bp identity matrix. The later setting can be viewed as 
adding subunit noise (Goldwynand Shea-Brown, 2011). Recently, 
it has been suggested that a different way of adding channel noise 
into the BSGs may result in more accurate stochastic behavior 
(Goldwyn and Shea-Brown, 2011; Goldwyn et al, 2011; Linaro 
et al, 2011; Orio and Soudry, 2012). The SDEs in (15) are of 
general form and do not preclude them. In fact, by setting B ! 
to be a block matrix with blocks equal to be the square root 
of the diffusion matrix for each ion channel, the channel SDE 
model (Goldwyn et al, 2011; Orio and Soudry, 2012) can easily 
be incorporated into (15). 

Finally, we note that, under appropriate technical conditions 
the SDE formulation applies to BSGs with voltage-gated ion chan- 
nels as well as other types of ion channels. The conditions require 
that the BSG model can be treated mathematically as a system of 
SDEs of the form (15) and that the latter satisfies the assumptions 
of Section 2.4.1. 

2.4.3. Overall encoding of the neural circuit model 

Taking into account the dendritic input from the feedforward 
DSPs and feedback DSPs, the encoding by the neural circuit 
model under the two noise sources is given by two systems of 
SDEs. With the Brownian motion Wj modeling the DSP white 
noise, the encoding of neuron i = 1,2, can be expressed as 



dY' = V (Y\ A) dt + B' (lA) dZ'(t), 



(16) 



where 



IZ 1 = 



v'dt + dW{ 
dWl 



dWl 



with v'(t) given by Equation (12). 



Note that in the system of Equations (16) the two output 
spikes trains (f^) , i = 1,2, k e Z, are the observables. Due to the 
intrinsic noise sources in the DSPs and in the BSGs, spike timing 
jitter may be observed from trial to trial by repeatedly applying 
the same stimulus to the neural circuit (see Section 2.3 in the 
Supplementary Material). 

3. ENCODING, DECODING, AND NOISE 

In Section 3.1 we present the mathematical encoding formalism 
underlying the neural circuit in Figure 1. We formulate stimulus 
decoding as a smoothing spline optimization problem and derive 
an algorithm that reconstructs the encoded signal in Section 3.2. 
Finally, we analyze the effect of noise on stimulus decoding in 
Section 3.3. 

3.1. ENCODING 

In this section, we formulate a rigorous stimulus encoding model 
based on the neural circuit shown in Figure 1. The input of 
the circuit is a signal U\ modeling a typical sensory stimulus as 
described in Section 2.2.1. The neural circuit generates a mul- 
tidimensional spike train that is assumed to be observable. We 
establish model equations by first describing the I/O relation- 
ship (i.e., the t-transform) of a single BSG. We then provide 
the t-transform of the entire neural circuit model that maps the 
input stimulus amplitude into a multidimensional spike timing 
sequence. 

3.1.1. The I/O of the BSG 

In the presence of a bias current P and absence of external 
inputs, each BSG in Figure 1 is assumed to be periodically spik- 
ing. Provided that the inputs are small enough, and by using the 
PRC, the BSG dynamics of spike generation can be described in 
an one-dimensional phase space (Lazar, 2010). 

Definition 3.1. A neuron whose spike times (f!), k e Z, i = 1, 2, 
verify the system of equations 

B' (x< (s - 4 + t' (s - 4, A) , Afj dZ\s) 

= r(r)-(A +1 -A), (m 

where 

dr> (f - 4, A) = [f (t - 4 + r' (f - 4, p) , p)f (18) 
B ! (x< (t - 4 + r< (t - 4, A) , r')) dZ'(t), 

with r' (0, P) = 0 and x' (f, P) the periodic solution to (13) with 
bias current P, is called a Project-Integrate-and-Fire (PIF) neuron 
with random thresholds. In (17), [•] denotes transpose and T'(P) 
is the period of limit cycle with bias current P. 
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As its name suggests, the PIF projects a weighted version of 
the input embedded in noise and the ion channel noise asso- 
ciated with the gating variables (B'dZ') onto the PRCs of the 
corresponding gating variables on a time interval between two 
consecutive spikes. Note that the integrand in (17) is identical to 
the RHS of (19). r'(f, V) on the LHS of (19) denotes the phase 
deviation and is driven by the perturbation on the RHS. The LHS 
of ( 1 7) represents the phase deviation measurement performed by 
the PIF neuron. The RHS of (17) provides the value of the mea- 
surement and is equal to the difference between the inter-spike 
interval and the period of the limit cycle. 

The BSG and the PIF neuron with random thresholds are, to 
the first order, I/O equivalent (Lazar, 2010). In Lazar (2010) it was 
also shown that a good approximation to the PIF neuron is the 
reduced PIF with random threshold. The functional description 
of the reduced PIF is obtained by setting the phase deviation in 
(17) to zero. 

Definition 3.2. The reduced PIF neuron with random threshold is 
given by the equations 



J2 f k+1 f n (s- i^^i Ufs-^Ay 

— 1 Jt 'k 

v')- -<£) +4 



(s)ds 



n = 1 1 
= T 



(19) 



where {el), k e Z, is a sequence of independent Gaussian random 
variables with zero mean and variance 



:[.;]» 

p = l Jt k n=\ 



ds. (20) 



For reasons of notational simplicity and without loss of general- 
ity, and unless otherwise stated, we shall assume here that B = I 
(N' = P'). The reduced PIF (rPIF) with random threshold can 
now be written as 



f +1 f[ ( s - 4, A v\s)ds = v (A - (4 +1 - 4) + 4, 

Jtl 



(21) 



where (e k ),k e Z, i= 1,2, is a sequence of independent 
Gaussian random variables with zero mean and variance 



E 



ds. 



(22) 



The above analysis assumes that the inputs are weak and therefore 
the BSGs operate on a limit cycle. Stronger signals can be taken 
into account by considering a manifold of PRCs associated with 
a wide range of limit cycles (Kim and Lazar, 2012). By estimating 
the limit cycle and hence its PRC using spike times, we have the 
following I/O relationship for each of the BSGs. 



Definition 3.3. The reduced PIF neuron with conditional PRC and 
random threshold is given by the system of equations 



ds ■ 



(23) 



where b' k = [V] 1 (t' k+1 — f£) , k e Z, is the total estimated bias 
current on the inter-spike interval [t' k , t' k+1 ]> I' 0 is an initial bias 
that brings the neuron close to the spiking region in the absence of 
input and (by abuse of notation) s' k ,k € Z, i = 1, 2, is a sequence 
of independent Gaussian random variables with zero mean and 
variance 



E 



ds, (24) 



and i/fj(s, b' k ) is the conditional PRC (Kim and Lazar, 2012). 

The conditional PRC formulation above allows us to sepa- 
rate BSG inputs into a constant bias current and fluctuations 
around it on short inter-spike time intervals. The bias current 
can be estimated between consecutive spikes, making the devi- 
ation from the limit cycle small in each inter-spike interval even 
for strong inputs. Moreover, by considering the conditional PRCs, 
the assumption that BSGs oscillate in the absence of input can be 
nearly dropped. Thus, it is not required for BSGs to always be 
on a limit cycle. Only when the neuron enters the limit cycle do 
we consider formulating the encoding using the rPIF model with 
conditional PRCs. 

Remark 3.4. Note that by parametrizing each of the PRCs with b' k , 
the variance of the error in (24) depends on the estimated PRC on 
each inter-spike interval. In conjunction with (23), we see that the 
variability of spike times depends on the strength of the input to the 
BSGs. 

3.1.2. The t-transform of the neural circuit 

The overall encoding by the neural circuit model can be 
expressed as 

f* +1 *[ ('-£*£) v '^ ds 

Substituting (12) into the above, we have 

f k+1 ir\ (s - t[, b[) [ h\ ll (s - r)u v (r)drds 

j t ' k+1 f[ (s - t' k , *4) jf h\ u (s -n,s- r 2 )uM) 



s' k , i = 1,2, k € Z. 



+ 



U\(r{)dr\dr2ds 
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Jt 'k 

- e r +i # ( s - £ ^ (^f) & - (25) 

l€Z Jt 'k 

-EE T +I ( s - 4. H) (^f ) (* - «/. * - f ™) * 
+ 4. «".;= i.2,« # j. 

We arrived at the following. 

Lemma 3.5. The model of encoding in Figure 1 is given in operator 
form by 



Ty k u\ + T 2 \u 2 = q[ + 4' ! ' = 1- 2, 6 Z, 



(26) 



w/iere mi g H}, «2 g 7i\, u 2 (ti, t 2 ) = u\(t\)u\(t 2 ), and, T^ k : 
7i\ — > M and 7^ : 7^ — > R «re bounded linear functionals 
given by 



Vk u i= f tk+] ( hi u (j-r)«i(r)<ifrfs, 
>/f[ JO] 

T 2 \ u 2 = / ¥>jfc(s) / ^2 h ( s_ r i' s ~ r 2)«2(n, r 2 )dr x dr 2 ds, 
Jt' k Jn 2 

4k = (H-^)/,%*w*-E r +1 ^( s )( p2/l f) 
EE (*-*/.* -4)*. 



^j(t-i^) 



(e[4] 2 ) : 



ande' k ,k e T,, are independent random variables with normal dis- 
tribution A/"(0, 1) and; = 1, 2, j ^ i. Equation (26) is caZZed f?ie 
t-transform (Lazar and Toth, 2004 ) of the neural circuit in Figure 1. 

Remark 3.6. X/ze t-transform describes the mapping of the input 
stimulus u\ into the spike timing sequence (tl), i=\,2,k G Z. 
r/jus, ("ne t-transform shows how the amplitude information of the 
input signal is related to or transformed into the time information 
contained in the sequence of output spikes generated by the neural 
circuit. 

We provide here further intuition behind the Equations (26). 
By the Riesz representation theorem (Berlinet and Thomas- 
Agnan, 2004), there exists functions <p\ k G H.\ such that 



(Ml, 



for all u\ G Tt\, 



and cf>' 2k G Ti\ such that 

T 2 \u 2 = («2, <P' 2 k)n\> for a11 "2 e ^2- 
Therefore, (26) can be rewritten in inner product form: 

0ijt> w i + («2> 02ife>Hi = Ik + s l (27) 

Recall that inner products are projections that are typically inter- 
preted as measurements. In the Equation (27) above, the signals 
u\ and u 2 are projected onto the sampling functions (j>U and <p' 2k , 
respectively. We also note that traditional amplitude sampling of 
a bandlimited signal U\ at times (f„), n G Z, can be expressed as 



(ui( • ), sinc(f„ - -)>l2( 



u\{t„), 



where sine (f) = sln ^ ^ is the impulse response of the ideal 
low pass filter bandlimited to Q} or in other words, the ker- 
nel of the RKHS of finite-energy band-limited functions (Lazar 
and Pnevmatikakis, 2009). Thus, the neural encoding model 
described by the Equation (27) can be interpreted as generalized 
sampling with noisy measurements with sampling functions <p\ k 
and <t>\ k . 

The formulation of the encoding model can easily be extended 
to the case when M neural circuits encode a stimulus in paral- 
lel. This is shown schematically in Supplementary Figure 1. A left 
superscript was added in the figure to each of the components to 
indicate the circuit number. 

3.2. DECODING 

In the previous section, we showed that the encoding of a signal u\ 
by the neural circuit model with feedforward and feedback DSPs 
and BSGs can be characterized by the set of t-transform Equations 
(26). We noticed that the Equations (26) are nonlinear in u\ due 
to the second order Volterra term. However, by reinterpreting the 
second order term as linear functionals on the higher dimen- 
sional tensor space TL\, (26) implies that the measurements taken 
by each of the neurons are the sum of linear measurements in two 
different vector spaces [see also Equations (27)]. 

In this section we investigate the decoding of signals encoded 
with the neural circuit in Figure 1. The purpose of decoding is 
to recover from the set of spike times the original signals, U\ (t) 
and «2 ( f l , h ) > that respectively belong to the two different vector 
spaces TL\ and Tl\. We formulate the decoding problem as the 
joint smoothing spline problem 



(u\, ii 2 ) = argmin 

«i eH\,u2<EH\ 



tl||«l 



+ A.2 1| «2 



i2 



+ EE( r iW+7>2-<) 2 



= Ijfc=l 



, (28) 



'Xkl-H 



where n' + 1 is the number of spikes generated by BSG i = 1,2. 
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Theorem 3.7. The solution to (28) is of the form 



M 2 (fl, fe) 



££44(ft.fe). 
>= i &= i 



where 
k= 1, 



,,(f) 



and 0i t (ti, *a) 



(29) 



1,2, 



Zi«ear equations 



4]' 



is ffoe solution of the system of 



( ($i + $ 2 ) 2 + + a 2 $ 2 ) c = (*i + * 2 ) q, 

<? 2 2 ] , flM d 



(30) 



$21 



$ 12 
<t> 22 



1, 2, 



and 



Proof: Proof of the theorem follows the Representer Theorem 
(Berlinet and Thomas-Agnan, 2004) and is given in detail in 
Appendix. 

Remark 3.8. When k\ = A 2 , the solution c amounts to 

c= ($! + 4> 2 + Aiir 1 q, 

where I is a« identity matrix of appropriate dimensions. 

Remark 3.9. Although (29) solves (28), in practice a minimum 
number of spikes is needed to obtain a meaningful estimate of the 
original signal. A minimum bound for the number of measure- 
ments/spikes can be derived in the noiseless case. Clearly, the bound 
has to be larger than the dimension of the space. This may require 
the signal to be encoded by a circuit with a larger number of neu- 
rons than the two shown in Figure 1 (Lazar and Slutskiy, in press). 
A number of such neural circuits in parallel can be used to encode 
input stimuli as shown in the Supplementary Figure 1. Theorem 3.7 
can be easily extended to solving the smoothing spline problem 



(«i, « 2 ) = argmin 

Ml €H{, u 2 eH\ 



Mll«lll^i + Ml «2 ll^i 



M 2 'V 



+ £££( mr >i + mT >2- m 4) 



m=li=li=l 



where m = 1,2, . . . , M, denotes the circuits number in 
Supplementary Figure 1. In addition, if the circuits consist of 
only first order feedforward kernels, then only U\(t) can be recon- 
structed. Similarly, if the circuits are comprised of only the second 
order feedforward kernels, then u%(t\, f 2 ) can be reconstructed but 
notui(t). 

Remark 3.10. Since « 2 (fi, f 2 ) = Mi(i"i)i/i(t 2 ) = w 2 (t 2 , ti), « 2 
belongs to a subspace of Til whose elements are symmetric func- 
tions. We also note that since the second order feedforward ker- 
nels are symmetric, the sampling functions (<p' 2k (t\, f 2 )) ,i= 1,2, 
k = 1, • • • , n', also belong to the same subspace. Therefore, if the 
sampling functions span the subspace of symmetric functions in H.\, 
U2 can readily be reconstructed with only (L l + l) (2L 1 + \) mea- 
surements/ spikes, rather than (2L 1 + l) 2 , the dimension oj r 7i\. 

Remark 3.11. The reconstruction of w 2 (fi,f 2 ) on D 2 strongly 
depends on the support (in practice the finite memory) of the kernels 
h\ u , i = 1,2 (see also Figure 2C). In the reconstruction example of 
the Supplementary Figure 5, we show that « 2 approximates m 2 well 
in the restricted domain where h\ u is nonzero. Outside this restricted 
domain, vanishes and « 2 is not well recovered as suggested by the 
large error in the Supplementary Figure 5E. 



3.3. EFFECT OF NOISE ON STIMULUS DECODING 

In this section, we investigate the effect of noise sources (i) on 
spike timing of the reduced PIF neuron, and (ii) on the decoding 
of stimuli encoded with a neural circuit. We will also present the 
effect of an alternative noise source model on both spike timing 
and stimulus decoding. 

3.3.1. Effect of noise on measurement and spike timing errors of the 
reduced PIF neuron 

As suggested by (22), the variance of the measurement error of the 
reduced PIF neuron is directly related to the PRC of the associated 
limit cycle. We first characterize the variance of the measure- 
ment error due to each individual noise source parametrized by 
the bias current I'. We then evaluate the spike timing variance 
between the spike trains generated by the Hodgkin-Huxley neu- 
ron and the reduced PIF neuron again as a function of the bias 
current P. We start with a brief description of the key elements of 
Hodgkin-Huxley neuron and the PIF neuron. 

We consider the stochastic Hodgkin-Huxley equations 



<«♦ = f (y\ A dt+dzxt), 



(31) 



where f is defined as in Section 2.2 of the Supplementary Material 
with additional normalization such that the unit of time is in sec- 
onds instead of milliseconds and the unit of voltage is in Volts 
instead of milivolts as conventionally used. Z'(f) takes the form 



dZ\t) 



l dt + o[dW[ 
a\dW[ 
a\dW\ 
a\dW\ 
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Here W' n (t) are independent standard Brownian motion pro- 
cesses and o' n ,n= 1,2,3,4, are associated scaling factors. 

The variance of the measurement error of the reduced PIF neu- 
ron due to each Brownian motion process W' n , n = 1, • • • , 4, is 
given by [see also Equation (22)] 

(e [4] 2 ) (?) = (rf [r„ (s - i, i 1 )] 2 ds. 02) 

We show in Figure 3A the variance of the measurement error 
in (32) associated with each source of noise of the reduced PIF 
neuron for the unitary noise levels a' n = l,n= 1,2, 3, 4. The 
variances given by (32) are plotted as a function of the bias current 

Clearly, the noise arising in dendritic stimulus processing ( W] ) 
induces the largest error, and together with noise in the potassium 
channels (W l 2 ), these errors are about two magnitudes larger in 
variance than those induced by the noise sources in the sodium 
channels {W' 3> W{). 

The above analysis is based on the analytical derivation of 
the measurement error in (32) for the rPIF neurons. The mea- 
surement error is closely related, however, to the spike timing 
variation of the BSGs subject to noise sources. A variance of 
1CP 6 in Figure 3A corresponds to a standard deviation of 1ms 
in spike timing. In practice the error between the spike times of 
the Hodgkin-Huxley neuron and the reduced PIF neuron can be 
directly evaluated. 



In order to do so, we randomly generated a weak bandlim- 
ited dendritic input. All evaluations were based on encoding a 
signal with the Hodgkin-Huxley neuron model described above 
with internal noise sources and bias current P. The spike times 
(t' k ) of the Hodgkin-Huxley neuron were recorded. Starting from 
each spike time tL we encoded the appropriate portion of the sig- 
nal by the reduced PIF neuron until a spike r t' k+l was generated. 
The difference between r t' k+l and t' k + 1 is the error in approximat- 
ing the encoding using the reduced PIF formulation. This process 
was repeated for each P. We computed the variance of the errors 
based on some 3000-5000 spikes generated in encoding the input. 

In Figure 3B, the variance of the spike timing error r t' k+1 — 
t' k+l for er n = 0, n = 1, 2, 3, 4, is shown. Since the reduced PIF 
is an approximation (even under noiseless conditions) and, 
although small, the error is nonzero. From Figure 3B, the vari- 
ance of the spike timing error is on the order of 10~ 9 . We shall 
evaluate the spike timing error variance of the intrinsic noise 
sources in a range much larger than 10~ 9 . 

We also tested to what extent each individual source of noise 
contributes to the variance of spike timing as suggested by the 
theoretical analysis depicted in Figure 3A. Indeed, the error vari- 
ance obtained through simulations in Figure 3C follows the basic 
pattern shown in Figure 3A. Figure 3C was obtained by setting 
one of the cr n 's to a nonzero value and the rest to 0 (the nonzero 
values were o\ = 02 = 0.01, er 3 = 04 = 0.1). Each nonzero value 
was picked to be large enough so that the error variance in the 




I 111 A/cm 2 ] I A/cm 2 ] 



FIGURE 3 I Variance of the measurement and spike timing errors. 

(A) Error measurement variances computed from the PRCs of the 
Hodgkin-Huxley neuron [Equation (32)]. Each individual variance is 
parametrized by the bias current /'. (B) Error variance between spike 
times generated by the noiseless Hodgkin-Huxley neuron and its 
reduced PIF counterpart. (C) The spike timing error variance due to 
each source of noise, obtained from simulations of the Hodgkin-Huxley 
neuron follow the pattern of the theoretically derived measurement 
error shown in (A). The spike timing error variances are obtained by 
setting, at each time, one of the <r n 's to a nonzero value and the rest 



to zero. The spikes generated by the Hodgkin-Huxley neuron are 
compared with the spikes generated by its reduced PIF counterpart. 
The variance of the differences between two spike times are 
normalized by the nonzero 07, mentioned before. (D) The spike timing 
variance due to the simultaneous presence of multiple noise sources 
approximates the sum of spike timing variances due to individual noise 
sources. Blue curve shows the spike timing variance obtained by 
simulating Hodgkin-Huxley equations using nonzero values for all 
oh, 1= 1,2,3,4. Red curve shows the sum of spike timing variances 
obtained in (C) with proper scaling. 
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absence of noise (Figure 3B) becomes negligible, and at the same 
time, it was small enough such that the states of the neurons did 
not substantially deviate from the limit cycle. To compare the with 
the ones in Figure 3A we normalized the error variance obtained 
in simulations by er„. 

Next, we tested whether the variance of spike timing due to 
presence of multiple noise sources is truly the summation of 
error variances due to individual noise sources. We simulated 
the Hodgkin-Huxley equations with o\ =02 = 0.005, 173 =174 = 
0.05. The total spike timing error variance shown in Figure 3D 
(blue curve) is very close to the sum of error variances in 
Figure 3C with proper scaling (red curve in Figure 3D). 

As suggested by the above analysis, the reduced PIF neuron 
with random thresholds largely captures the encoding of stimuli 
by BSGs subject to intrinsic noise sources. 

3.3.2. Effect of noise on stimulus decoding 

In order to quantitatively explore how noise impacts signal 
decoding, we recovered from spikes the signal encoded by the 
noisy neural circuit of Supplementary Figure 1. We started with 
the base-level noise-less case described in Section 3.2 of the 
Supplementary Material (M = 4) and proceeded to introduce 
individual noise terms with a range of scaling factors. For exam- 
ple, we set a l 2 = = a\ = 0 and varied a\. We also tested the 
case when 10<7j = IOct^ = o\ = o\ for the aggregated effect on 
stimulus recovery. We choose to use a\ and a\ 10 times larger 
than a\ and a\ so that each noise source introduced a similar 
error. 

In all simulations, the Euler-Maruyama scheme (Kloeden and 
Platen, 1992) was used for the numerical integration of the 
SDEs. We performed 20 encoding and decoding experiments. 
Each time, the input stimulus was generated by randomly pick- 
ing from a Gaussian distribution the real and imaginary parts of 
the coefficients u; in (1). We further constrained the stimuli to be 
real-valued. (An example is given in Supplementary Figure 5.) For 
each noise level, the input signal was encoded/decoded. The mean 
Signal-to-Noise Ratio (SNR) across 20 experiments is reported 
for each noise level. The SNR for the reconstruction of u\ was 
computed as 



SNR = 10 log! 



INI 



\U\ — U\ I 



(33) 



where u\ is the original signal and u\ is its reconstruction. Note 
that the spike time occurrences generated for the same signal are 
different for each noise level. Since the sampling functions are 
spike time dependent, the number of measurements/spikes may 
not be the same for each noise level. Moreover, at times, the sam- 
pling functions may not fully span the stimulus space. To reduce 
the uncertainty caused by the stimulus dependent sampling we 
averaged our SNR data over 20 different signals. 

Figure 4A shows the SNR of the reconstruction of signal 
U\{t) against different noise strength. Figure 4B shows the SNR 
of the reconstruction of signal u\{t) = ui{t, t). The reconstruc- 
tion SNR in Figure 4A largely matches the inverse ordering of 
noise strength of each of the individual noise sources shown 
in Figure 3A. DSP noise sources degrade the reconstruction 



performance most strongly while noise sources associated with 
gating variables m and h have a much smaller effect for the same 
variance level. Since the variance of measurement error is the sum 
of error variance in each variable, the case when lOeri = 10(72 = 
a 3 = a 4 = a exhibits the lowest performance. 

3.3.3. Effect of an alternative noise model on spike timing and 
stimulus decoding 

Biologically, the effect of channel noise on the operation of the 
BSGs is due to the ON-OFF activity of a finite number of ion 
channels. The Hodgkin-Huxley equations and the noise terms 
used in Section 3.3.2 correctly capture the dynamics in the limit 
of infinitely many channels. Recent research, however, suggests 
that the model equations may not correctly model the ion cur- 
rent fluctuations for a finite number of channels (Goldwyn and 
Shea-Brown, 2011). 

We consider here an alternative stochastic formulation of the 
Hodgkin-Huxley model that more precisely captures the ion 
channel kinetics. By using a finite number of ion channels the 
strength of noise amplitude becomes directly related to the actual 
number of ion channels. Therefore, the free variables are only the 
number of potassium and sodium channels that are both biologi- 
cally meaningful. The successful use of an alternative noise model 
as described below also suggests that our analysis can be applied 
to a wide range of stochastic formulations of BSGs based on SDEs. 

We shall construct here stochastic ion channels using conduc- 
tance noise rather than subunit noise as in the previous Sections 
(Goldwyn and Shea-Brown, 2011; Goldwyn et al, 2011). This 
stochastic Hodgkin-Huxley system is simulated using a diffusion 
approximation following (Orio and Soudry, 2012). The system of 
SDEs can be expressed as 

dY' = f (j 1 , 1 1 ) dt + B l (Y>) dZ'(t), 

where Y 1 has 14 state variables and the full system can be found 
in Section 3.3 of the Supplementary Material. Here i = 1 for 
simplicity. 

The variance of the measurement error is now given by (20). 
We can decompose the variance into three terms as 



E 



[4] 2 = e[4 v ] 2 + e[4J 2 + e[4 N(! ] 2 



where e' kv , s' kK , e' kNa are measurement errors associated with the 
noise in the DSP, in potassium channels and in sodium channels, 
respectively. 

As e' kv is quantitatively the same as that in Section 3.3.2, we 
focus here on e' kK and e\ Na - The variance of the errors can be 
respectively expressed as 



;[4j J to 

e r +i fe r - {s - **■ nb "p ( x,(s - **■ n ) 

, T ^1. 1 



ds. 
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FIGURE 4 | SNR reconstruction error of encoded signals with a total of 
M = 2 circuits (4 neurons). Color legend: (Blue) a\ = a, a' 2 = ai = = 0. 

(Green) a' 2 = a, a\ = = = 0. (Red) 0-3 = a, ai = = = 0. (Black) 
= a, g -j = = CTj 1 = 0. (Magenta) IOct] = 10cr^ = trj = 04 = a . In-sets (on 



the left) are typical reconstructions that yield corresponding SNR indicated by 
arrows. The top left in (A) shows an example of reconstruction (green) whose 
SNR is 25 dB when compared to the original signal (blue). (A) SNR of 
reconstruction of Ui (t). (B) SNR of reconstruction of u\(t) = U2(t, t). 



and 



11 



ds. 



Note that b n p, n = 1, • • • , 14, p = 2, 3, • • • , 15, are functions 
that dependent on either the number of potassium channels N^a 
or the number of sodium channels Njc, and the states of the 
neuron. 



We first evaluate (E [e' kNa ] )(/') using the PRCs. The PRCs 
are obtained by letting Nxa = Nk = 00 and thereby making the 
system deterministic. Since the measurement error variance for 
fixed I' is proportional to (N^ fl ) _1 , it is shown in Figure 5A 
as a function of the bias current /' for Nn„ = 1. Similarly, the 

variance of the measurement error [e^] 2 ^ (P) for Nk = 1 is 

shown in Figure 5 A, and it is proportional to (Afc) -1 for a fixed 
P. We notice that, when the number of channels is the same, 
the measurement error due to the sodium channels is of the 
same order of magnitude with the measurement error due to the 
potassium channels. However, the number of sodium channels is 
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FIGURE 5 | The variance of the measurement and spike timing error 
associated with the sodium channels (blue) and the potassium 
channels (red) of the Hodgkin-Huxley equations with alternative noise 
sources parametrized by the bias current /. (A) The variance of the 
measurement error computed from PRCs of Hodgkin-Huxley equations, 
with N Ng = 1 and N K = 1 . Actual variance with different number of ion 
channels is inversely proportional to Nn a and Nk, respectively. (B) Spike 
timing variance obtained in simulations by comparing the spike times 
generated by the Hodgkin-Huxley neuron with channel noise and the spike 
times generated by its reduced PIF counterpart. Blue curve is obtained by 
using N Na = 5 x 10 4 , N K = oo, and normalized to 1 sodium channel. Red 
curve is obtained by using Nk = 5 x 10 4 , N/v a = oo, and normalized to 1 
potassium channel. 



typically 3-4 times larger than the number of potassium channels. 
Therefore, in contrast to the previous case, the error induced by 
sodium channels shall be larger than that induced by potassium 
channels. 

We also analyzed in simulations the difference between spike 
times generated by the alternative stochastic formulation of the 
Hodgkin-Huxley equations and those generated by the corre- 
sponding reduced PIF neuron. We used in simulation ATjva = 

5 x 10 4 , Nk = oo, to obtain the variance (k [e* kNa \ 2 ^ (P) and 
scaled it by IVjv fl to compare it with Figure 5A. Similarly, we used 
N K = 5 x 10 4 , N Na = oo, to obtain the variance E [sj^] (P). 
The spike timing variances of error across different I' are shown 
in Figure 5B The pattern of similarity between variances in 
Figures 5A,B suggest that the reduced PIF with random threshold 
associated with this formulation of Hodgkin-Huxley equations is 
highly effective in capturing the encoding under internal noise 
sources. 

We now show how ion channel noise sources affect the decod- 
ing of the input signal. We varied the number of sodium chan- 
nels Nmu and fixed the number of potassium channels to be 
Nk = 0.3NjVa, a ratio typically used for Hodgkin-Huxley neurons 
with the alternative noise source model. By decoding the input 
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FIGURE 6 | SNR of reconstruction of the signals. (A) SNR of u, (f). (B) 
SNR of uj(r) = u 2 (r, f). (Blue) N Nb = N, % = 0.3W Ns . (Green) N Ng = N, 
N K = oo. (Red) N Nb = oo, N K = N. 



stimulus we show how increasing the number of ion channels 
improves the faithfulness of signal representation. The SNR of the 
reconstruction of U\(t) and u\{t) are depicted in Figure 6. SNR 
goes down to about 4dB when 600 sodium channels and 180 
potassium channels are used. This corresponds to a membrane 
area of about 10 \im 2 with a density of 60 |xm 2 in sodium chan- 
nels and 18 (im 2 in potassium channels (Goldwyn et al., 2011). 
We also tested the reconstruction for the case when one type of 
ion channels is infinitely large, i.e., deterministic, while varying 
the number of ion channels of the other type. The result is also 
shown in Figure 6. The noise from the dendritic tree shall have 
similar effect on the representation since the voltage equation is 
the same as in Section 3.3.2. 

4. FUNCTIONAL IDENTIFICATION AND NOISE 

In Section 4.1 we show how to functionally identify the feedfor- 
ward and feedback DSPs of the circuit in Figure 1 under noisy 
conditions. We assume here that the BSGs have already been iden- 
tified using a methodology such as the one developed in Lazar 
and Slutskiy (2014). In Section 4.2 we discuss the effect of noise 
parameters on the quality of DSP identification. 

4.1. FUNCTIONAL IDENTIFICATION 

In the decoding problem analyzed in Section 3.2, we recon- 
structed unknown input stimuli by assuming that the neural 
circuit in Figure 1 is known and the spike trains are observable. 
In contrast, the objective of the functional identification prob- 
lem investigated in this section is to estimate the unknown circuit 
parameters of the feedforward and feedback DSPs from I/O data. 
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The I/O data is obtained by stimulating the circuit with control- 
lable inputs and by measuring the time occurrences of the output 
spikes. This basic methodology has been a standard practice in 
neurophysiology for inferring the function of sensory systems 
(Hubel and Wiesel, 1962). We assume here that either the BSGs 
are known in functional form or the family of PRCs associated 
with the BSGs have already been identified (Lazar and Slutskiy, 
2014). 

Although decoding and functional identification are seemingly 
two different problems, they are closely related. By exploiting the 
commutative property of linear operators, we can rearrange the 
diagram of the neural circuit model of Figure 1 into the form 
shown in Figure 7. We notice that the outputs of Figure 7 and 
those of Figure 1 are spike time equivalent, as long as the RKs 
K 2 and Ki have large enough bandwidth. In what follows we will 
evaluate the four Volterra terms, i.e., the four dendritic currents 
feeding the BSG of Neuron 1 in Figure 7. 

Formally, the first order (feedforward) Volterra term can be 
written as (Lazar and Slutskiy, in press) 

f h\ 1, (t-s)u l (s)ds= f M i(f-s)(P 1 /j} 1 ')(s)ds. (34) 

Similarly, the second order (feedforward) Volterra term 
amounts to 

/ h\ h (t - si, t - 52) u 2 (si, 52) ds\ds 2 (35) 

Jd 2 



= f u 2 (t-s 1 ,t-s 2 )(v 1 h l 2 1, )(s l ,S2)ds 1 ds2. 

JV>2 V ' 

The above equations suggest that the projections of the feedfor- 
ward kernels can be re-interpreted as inputs, whereas the signals 
u\ and «2 can be treated as feedforward kernels. 

In Section 2.2.2 we introduced two RKHSs, Ti. 2 and 7i\, for 
modeling two different spaces of spikes. The elements of TL\ are 
functions defined over the domain [0, S 2 ] with 

S 2 >su PP )^) + max)( f ; +1 -^))^ i2teZ . 

The period S 2 is large enough to ensure that any spike that arrives 
supp{/ij ; '} seconds prior to the arrival of tl, or earlier, will not 
affect the output of the feedback kernel on the inter-spike time 
interval [t' k , t' k+1 ]. Thus, such spikes will not introduce additional 
error terms in the f-transform evaluated on the inter-spike time 
interval [t' k , fjj. , J. Note that the domain [0, S 2 ] of the functions 
in Ti. 2 may not be the same as the domain of the input space Ji \ . 
However, such a domain can be shifted on a spike by spike basis to 
[fjj , l — S 2 , t' k+ J for the inter-spike time interval [t' k , fjj , J. This 
is important for mitigating the practical limitation of modeling 
the stimuli as periodic functions in 7i\ . 

The response of the first-order feedback term to its spik- 
ing input on the inter-spike time interval [fjj., t' k+1 ] in Figure 7 
amounts to (i / j) 
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FIGURE 7 | Diagram of the neural circuit that is spike timing 
equivalent with the one in Figure 1 highlighting the duality 
between neural decoding and functional identification. Note that 
the input stimuli and the DSP projections are reordered to reflect 



that the unknowns are the DSP projections. The input stimuli u](t), 
U2(fi,f2), and the kernel representation of spikes (see also 
Section 2.2.2) are intrinsic to the neural circuit. The DSP projections 
are interpreted as inputs. 
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(v 2 h 2 fj (t) * ^2 0 = X! (T 1 ^)^ - t|) for identification, the t-transform (26) can be written as 



It is clear from Section 2.2.2 that 



E ( p2/! f) (t 



4> 



(36) ^"[P'/z} 1 '] + m C 1 2 i k [V 1 h 1 2 li ] + m £f k [V 2 hf] + m C 2i k [P 2 hf] 
= m <t k + m ei (38) 

hfit-fy for i,j=i,2,i jfe;,jk € Z. Here m £jj. : H\ -* R, : 



(4 +1 -s 2 <<W-, 



if is at least larger than the effective bandwidth of \xf and 
L 2 -> oo. 

Similarly, the response of the second-order feedback kernel 
to its spiking input on the inter-spike time interval [t' k , t' k+l ] 
amounts to 

E E (v 2 hf)(t-it-ti) 

ttl + 1 -S i <t><tl + i n:t- c+l -!?<ti,<tl +l 



E 



^ +1 -s 2 <^<4 +1 »4 +1 -s 2 <^<4 +1 



J2 hf(t-it-i) 



(37) 



if £2 2 is large enough and L 2 -> oo. 

Combining (34), (36), (36), and (37), for each spike interval 
[t' k , t' k+ J, the aggregated output current of the DSPs of Neuron z 
in Figure 7, shall converge to the DSP aggregated output current 
of Neuron i in Figure 1 for large enough Q, 2 . A direct consequence 
of this equivalence is that, under the same additive Gaussian white 
noise and channel noise in the BSGs, the t-transform of the circuit 
in Figure 7 and in Figure 1 are identical. 

Note that the outputs of the feedforward kernels are always 
equivalent; the equivalence of the outputs of the feedback ker- 
nels requires, however, the use of large enough bandwidth £1 2 . 
Otherwise, the equivalence in the t-transform is invalid and an 
additional noise term appears in the t-transform of the Neuron 1 
in Figure 7. 

The projections of the Volterra DSP kernels of Figure 7 are 
interpreted as inputs, while the input stimuli and the train of 
RKs at spike times replace the impulse response of the corre- 
sponding filters. Therefore, the functional identification problem 
has been transformed into a dual decoding problem, where the 
objects to decode are the set of projections of Volterra DSP ker- 
nels and the neural circuit is comprised of "stimulus DSP kernels" 
and "spike DSP kernels" with the same BSGs and noise sources. 
The only difference is that, instead of a Single-Input Multi- 
Output decoding problem, the identification was transformed 
into a Multi-Input Multi-Output decoding problem. In addition, 
multiple trials using different stimuli are needed; this procedure 
is illustrated in block diagram form in Figure 8. By stimulating 
the neural circuit with multiple stimuli in the functional iden- 
tification setting, multiple neural circuits effectively encode the 
projections of the DSP kernels. 

We are now in the position to derive the t-transform of Neuron 
1 in Figure 7. Assuming that m = 1 , • • • , M, trials are performed 



7i\ R are bounded linear functional associated with the 
feedforward DSP kernels, and m Cf k :H 2 ^R, m Cf k :H 2 ^-R 
are bounded linear functionals associated with the feedback DSP 
kernels for each trial m. The above functionals are defined as 



m jC] i k [V 1 h\ li ] = 



£ 2k [T h 2 ] 



C lk YP h 1 ] 



n m t l P 

I k + lm 4(s) I u\'{s-r)(V l h\ li ){r)drds 
J m f k Jni 

/ k+lm 4(s) / « 2 m (5- 

J m t{ Jo 2 



r u s-r 2 ){V l h\ u ) 



(n, ri)dr\dnds, 



E / , 



<p' k (s)(V 2 hf) 



L T k+1 a - t l < v k+l 



!-S 2 < "«{ < , "■■'" t 'k+ I" 5 ' =£ m 'n < m A+l 1 

[">'(*) • (V 2 hf)( S -%S- m i)]ds 
q< = CX -<*/<•) / k + lm V [{s)ds, 
ir[{t- m i, m V k ) 



>i(f) = 



and m e' k , i=\,2,k € Z, m = 1, • • • , M, are independent ran- 
dom variables with normal distribution 7V(0, 1). 

The functional identification of the neural circuit in Figure 7 
can then be similarly defined to the decoding problem. We for- 
mulate the identification of the noisy neural circuit again as two 
smoothing spline problems, one for each neuron, 



(V l h\ n ) 

(p^f 1 ) 



argmin 



lp=l,2 



£X>?lrt Pl || 2 (39) 
p=l r=l 



M "V 111 \ 
m=ljfc=l U = lr=l J 
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FIGURE 8 | Diagram of the functional identification with multiple trials. The neural circuit is presented a different stimulus u™(f) for each trial m. See also 
Figure 7 for details of a single trial. 
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and 



(T^hf 2 ) 
(T^hf 2 ). 



argmin 



\VW en 



J/>=1,2 



J^E^hr W 2 C40) 

p= 1 r= 1 



M "V / 2 2 



+ 



EE EE m ^ u ] 



m=lk=l \p=l r=l 



where '"n, is the number of spikes generated by Neuron i in 
trial m. 

The solution can be obtained in a similar way as in 
Theorem 3.7. 

Theorem 4.1. The solutions to (40) is 

M '"n 1 



m = 1 k = 1 

M '"n 1 

(^"Xfi, fe) = E E m ^ m fe(fi, h) 

m=\k=\ 

M '"n 1 

Cp2fc221 )(f) = ^ ^ ™ c ,"'0 3 ,(f) 

m=\k=\ 

M "'n 1 

(r 2 hi 21 xtu f 2 ) = E E m <* m &*(i. fe), 

m= 1 t= 1 



where 



c = ['ci • • • 1 ci „i , 



, • • • , ci • • • CM„lJ , 
is the solution to the system of linear equations 

((*! + $2 + *3 + *4) 2 + + X\<S> 2 + 4*3 + A|* 4 ) C 



(*1 + $ 2 + *3 + *4)q, 



(41) 



where 



and 



Hi HM„i] 



1M -i 



$ M1 . . . $ MM 



[*r] H = < m *a."*a>. 

In addition, the sampling functions m (j>ik are given by 
(pxkit) = L. lk K l]v 
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Proof: The proof is similar to the one of Theorem 3.7. □ 

Since each of the kernel projections may be in a different 
RKHS, and their domain may also be different, the identification 
of all filters resemble that of the multi-sensory Time Encoding 
Machines. Recall that multi-sensory TEMs encode within the 
same circuit time-varying and space-time varying sensory signals 
while decoding remains tractable (Lazar and Slutskiy, 2013). The 
solution to (41) can similarly be obtained as the solution to (40) 
above. 

Note that we are only able to identify the projection of the 
Volterra kernels. This is because, in practice, we can only probe 
the system with signals in a bandlimited space. By increasing the 
bandwidth of the elements of the Hilbert space, the projection 
of the kernels will converge to their original form (Lazar and 
Slutskiy, 2012). 

Remark 4.2. It is important to note that in order to have a good esti- 
mate of the kernels, stimuli must fully explore all input spaces. This 
can he quite easily achieved for the feedforward DSP kernels by using 
many (randomly generated) signals that cover the entire frequency 
spectrum. However, to properly identify the feedback DSP kernels, 
spike trains must be diverse enough to sample its different frequency 
components. This may not be easy to realize in practice. For first 
order feedback kernels, spike trains with constant spike intervals are, 
for example, undesirable. The Fourier transform of regular Dirac- 
delta pulses is a train of Dirac-delta pulses in the Fourier domain. 
This means that only certain frequency responses of the DSP kernels 
are, for example the DC component, sampled. The rest of the fre- 
quency components are essentially in the null space of the sampling 
functions m (f>ik, i= 1,2, m = 1, • • • , M. Similar effect applies to 
the space of trigonometric polynomials. If the spike intervals exhibit 
small variations, many of the frequency components may be sam- 
pled but the energy at DC may be too dominant. In this case, noise 
may contaminate more severely the measurement of non-DC com- 
ponents and may yield unsatisfactory identification. This effect may 
pose even more stringent requirements on the identification of the 
second order feedback kernels, as it requires the interaction between 
two spike trains. 

42. EFFECT OF NOISE ON IDENTIFICATION 

In order to evaluate the effect of noise on the identification of 
the neural circuit in Figure 1 we included intrinsic noise into 
the example neural circuit discussed under noiseless conditions 
in Section 4.1 of the Supplementary Material. Randomly gener- 
ated signals were used in the identification examples given here. 
Chosen in the same way as in the decoding example in Section 
3.3.2 all these signals are used here to identify the neuron in ques- 
tion. Therefore, in this section, multiple signals are used in repeat 
experiments to identify the parameters of a neural circuit. By 
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contrast in Section 3.3.2, multiple neurons are used to encode a 
single signal. 

First, we evaluated the effect of noise on the quality of iden- 
tification of DSP kernels of Neuron 1 in Figure 7 with a BSG 
modeled by the SDE (31) with IOctJ = 10<rj = a\ = a\ = a. 
Figure 9 shows the SNR of the identified DSP kernels in Figure 7 
across several noise levels a. As expected, the general trend for 
all four kernels is decreasing SNR with increasing noise levels. We 
notice that the identified feedforward DSP kernels have similar 
shape as the original kernel, even at high noise levels. However, 
the feedback DSP kernels undergo a change in shape at high noise 
levels. We can see that the time instance of the peak amplitude 
in the first order feedback kernel is shifted to an earlier time 
instance. 



Second, we investigated the identification of DSPs for a BSG 
noise model already described in Section 3.3.3. Figure 10 shows 
the SNR of the identified DSP kernels across a different number of 
sodium channels N^ a while Nk = 0.3Nmu- The SNR plots suggest 
that the identification quality increases as more ion channels are 
present in the BSGs. 

Additionally, as discussed in Remark 4.2, BSG noise sources 
may degrade severely the identification of feedback kernels when 
the spike trains generated in trials are not sufficient for explor- 
ing the two spike input spaces. We show an example of the later 
in Figure 11. The two BSGs have higher bias currents and lower 
input current magnitude. The later was achieved by scaling down 
the magnitude of the DSP kernels. The combined effect results 
in regular spiking intervals in both neurons. The identification 




FIGURE 9 | SNR of identified DSP kernels. Noise added using SDE (31), with 
IOctJ = IOctj — (To = tri — <r. (A) Kernel /i] 11 . In-sets provide a comparison 
between the original and the identified kernel. (B) Kernel h" 1 . In-sets are 



identified kernels. Original kernel is on the lower left. (C) Kernel ft 221 . In-sets 
provide a comparison between the original and the identified kernel. (D) 
Kernel h?, 21 . In-sets are identified kernels. Original kernel is on the lower left. 
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result under noiseless conditions is shown in Figure 11. Note that 
since the t-transform of the Hodgkin-Huxley neuron is not exact, 
a small error is introduced even if intrinsic noise is not present. 
We see that the feedforward DSP kernels can be identified quite 
well, yielding SNRs of around 17 dB. However, the feedback DSP 
kernels are not well identified. In particular, the identified second- 
order feedback kernel has a wide spread, similar to the high noise 
case in Figure 9D. This suggest that the spike pattern is not suffi- 
ciently exploring the space of feedback kernels. A large number of 
frequency components are only weakly sampled and they can be 
very easily contaminated by noise. The presence of both intrinsic 
noise sources can exacerbate the condition further. This effect is 
confirmed with a simulation of the circuit using Integrate-and- 
Fire (IAF) neurons. Since the f-transform for the IAF neuron is 
exact (Lazar and Toth, 2004), both feedback kernels can be iden- 
tified even if the generated spikes only weakly explore certain 
frequency components. However, by artificially adding a small 
measurement error to the t-transform of the circuit with IAF neu- 
rons, similar results to those in Figure 11 can be obtained (data 
not shown). 

5. DISCUSSION 

In this paper, we introduced a novel neural circuit architec- 
ture based on a neuron model with a biophysical mechanism 
of spike generation and feedforward as well as feedback den- 
dritic stimulus processors with intrinsic noise sources. Under 
this architectural framework, we quantitatively studied the effect 
of intrinsic noise on dendritic stimulus processing and on 
spike generation. We investigated how intrinsic noise sources 
affect the stimulus representation by decoding encoded stim- 
uli from spikes, and quantified the effect of noise on the 
functional identification of neural circuits. We argued that a 
duality between stimulus decoding and functional identifica- 
tion holds. Therefore, the encoding framework based on the 
neural circuit architecture studied here can be applied to both 



the reconstruction of the encoded signal and the identifica- 
tion of the dendritic stimulus processors. We systematically 
showed how the precision in decoding is affected by differ- 
ent levels of stochastic variability within the circuit. These 
results apply verbatim to the functional identification of den- 
dritic stimulus processors via the key duality property mentioned 
above. 

Our theoretical framework highlights two indispensable com- 
ponents of modeling signal representation/processing in a neural 
circuit — dendritic stimulus processing and spike generation. Such 
a divide and conquer strategy is ubiquitous in engineering circuits 
and leads to a separation of concerns. Recent experimental stud- 
ies also showed that interesting nonlinear processing effects take 
place first in the dendritic trees rather than in the axon hillock 
(Yonehara et al, 2013). 

We presented here two types of nonlinear dendritic stimulus 
processors. The first type are feedforward DSPs that respond to 
continuous input sensory stimuli. The second type are feedback 
DSPs that respond to spiking inputs. We quantitatively demon- 
strated how intrinsic noise sources would affect the identification 
quality of all these DSPs. The examples in Section 4.2 suggest that 
in identification feedback kernels are more vulnerable to internal 
noise sources than feedforward kernels. In particular, the overall 
shape of the identified feedback kernels differs significantly from 
that of the underlying kernels when the strength of noise sources 
becomes large. Meanwhile the identified feedforward kernels are 
qualitatively preserved, albeit not accurately. 

Most of the single neuron models (LIF, QIF) in the literature 
have focused on the spike generation mechanism. The encoding 
capability of these models is typically investigated based on rate 
encoding (Eliasmith and Anderson, 2003; Lundstrom et al., 2008; 
Ostojic and Brunei, 2011). For both decoding and identification 
we used here the occurrence times of spikes generated by spik- 
ing neuron models. Most importantly, the BSG models discussed 
here were characterized by a PRC manifold (Kim and Lazar, 2012) 
in the presence of noise, while many simplified models (such as 
the LIF) can be effectively described with a single PRC. Other 
sensory neuron models, e.g., GLM (Pillow et al, 2011), usually 
rely on a rate-based output or Poisson spike generation that do 
not take into account key advances in dynamical systems-based 
spiking neuron models. 

As already mentioned before, we investigated how intrin- 
sic noise sources affect the stimulus representation by decoding 
encoded stimuli from spikes. We are not suggesting, however, 
that the decoding algorithm considered here is implemented in 
the brain. Rather, we argue that decoding is effective in measur- 
ing how well information is preserved in the spike domain. The 
decoding formalism allowed us to investigate how noise affects 
the fidelity of signal representation by a population of neurons by 
reconstructing stimuli and comparing their information content 
in the stimulus space. 

While decoding can serve as an "oscilloscope" in understand- 
ing stimulus representation in sensory systems, functional identi- 
fication serves as a guide in experiments to functionally identify 
sensory processing. Based on spike times, the identification algo- 
rithm presents a clear bound on the number of spikes that are 
necessary for perfect identification under noiseless conditions. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



September 2014 [ Volume 8 | Article 95 | 20 



Lazar and Zhou 



Neural circuits with intrinsic noise sources 



20- 



= 




-20- 



— original kernel 

— identified projection 



0.1 



0.2 
time [sec | 



0.3 



0.4 



20 



-20 



1 



— original kernel 

— identified projection 



0.05 0.1 
time |sec] 



0.15 



0.02 0.04 0.06 0.08 
tj [sec] 



0 0.02 0.04 0.06 0.08 
tj [sec] 



0 0.02 0.04 0.06 0.08 
tj Isec] 



100 




FIGURE 11 | Examples of functional identification when the 
generated spikes do not fully explore the space of feedback 
kernels. (A) Original first order feedforward kernel (black) and identified 
projection of the kernel (red). (B) Original first order feedback kernel 
(black) and identified projection of the kernel (red). (C) Original second 



order feedforward kernel. (D) Identified projection of second order 
feedforward kernel. (E) Error of identified second order feedforward 
kernel. (F) Original second order feedback kernel. (G) Identified 
projection of second order feedback kernel. (H) Error of identified 
second order feedback kernel. 



Phrased differently, when a certain number of spikes are acquired 
from a neuron of interest, the identification algorithm places 
a constraint on the maximum DSP kernel bandwidth that can 
perfectly be recovered. 

In more practical terms, we advanced two important applica- 
tions of the circuit architecture considered in this paper. The first 
one considers dendritic stimulus processors that process infor- 
mation akin to complex cells in VI. The second one adapts the 
widely used Hodgkin-Huxley model known in the context of 
neural excitability (Izhikevich, 2007) and analysis of neuronal 
stochastic variability to stimulus encoding in the presence of 
noise. 

Based on the rigorous formalism of TEMs (Lazar and Toth, 
2004), we extended our previous theoretical framework (Lazar 
et al., 2010) and argued that spike timing is merely a form 
of generalized sampling of stimuli. By studying sampling (or 
measurements) in the presence of intrinsic noise sources, we 
showed to what extent neurons can represent sensory stimuli in 



noisy environments as well as how much noise the identification 
process can tolerate while preserving an accurate understanding 
of circuit dynamics. 

The reconstruction and identification quality are certainly not 
only related to the strength of noise, but also the strength of the 
signal. In particular, when the signal strength is small, two facts 
may affect the quality of reconstruction. First, neurons may not 
produce enough spikes that have valid f-transforms. Second, they 
may be contaminated by even weak noise, i.e., the signal-to-noise 
ratio is low. It is well known, however, that neural systems use gain 
control to boost the relevant signal (Shapley and Victor, 1978; 
Wark et al., 2007; Friederich et al, 2013). Such strategy may be 
useful for increasing the signal strength relatively to the strength 
of the noise. Gain control may also suppress large signals to fit 
into the range of operation of the BSGs. The gain control itself, 
maybe considered as a type of Volterra feedforward DSP kernel 
(Lazar and Slutskiy, in press) and the interaction with feedback 
loops driven by spikes. The lack of spikes may be compensated 
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by adding other neurons that are sensitive to other features in the 
input stimuli. 

A key feature in our neural circuit model is the nonlinear 
processing in the feedforward and feedback paths. Nonlinear 
interaction between feedforward DSPs and feedback DSPs have 
not been considered here. However, they are of interest and could 
be addressed in the future. Self-feedback was not included in 
the model for clarity, but can be considered within the frame- 
work of our approach. Self-feedback was introduced to add 
refractoriness to phenomenological neuron models (Keat et al., 
2001; Pillow et al, 2008). Our BSC models, on the contrary, are 
conductance-based models that have a refractory period built in. 

Throughout this paper we assumed that the BSGs themselves 
have been perfectly identified. The intrinsic noise in the BSGs may 
degrade the identification quality of conditional PRCs. This may 
result in a lower identification quality as shown in the examples. It 
is beneficial to investigate in the future a method that can identify 
the entire circuit at once so that intrinsic noise in the circuit only 
affects the identification process a single time. 

The theoretical results obtained here suggest a number of 
experiments in the early olfactory system of fruit flies. The 
glomeruli of the antennal lobe can be modeled using the Volterra 
DSPs discussed here and the projection neurons in the anten- 
nal lobe are accessible by patch clamping (Lazar and Yeh, 2014). 
Functional identification of DSPs can then be carried out for 
studying olfactory stimulus processing in an accessible circuit 
with intrinsic noise sources (Masse et al, 2009). 
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that 73, «i = (Ml > 0u) i V«i € 7Y}. Moreover by the reproducing 
property 



"20' 

orthogonal, with 



are 



H\ 0 = span ([4>i k }, i = 1,2, fc = 1, • • • , rA 

4>\ k {t) = [4>[ k ,K\}j=T^ t . and 



H\k = f«2 e H|73 t « 2 = 0, i = 1, 2, k = 1, 
Finally, 



Let TL } 0 be a linear subspace of ?ij spanned by 0^. 

Wio = span (j^}, i = 1,2, k= 1, •••,«') 2 ^ 

,, X XI (-^"lO + ^2* M 20 _ Ik) + MIN0IIL1 + ^-2 II "20 II 

and let 7ijg be a linear subspace of TL[ defined by j=lk=1 1 

H\in = lui e H\\Tj k ui = 0,i= l,2,k= l,--- V^Y^iVw,, ,„1utu„ _i_„-^ vV 



= lit=l 



Then, for any Mi 6 H\£ and any £ 2 = , £jf = i 4^1* e W !o' we + Xi||m 10 ||^ , + X 2 ||m 2 oI| 2 i 

have "i w 2 

■ E E 4rf*) = E E = E E ««i = «• 5 £ S ( ^ (U1 ° + 4) + ^ (M20 + 4) " ^ 

j=lt=l ' i = 1 fc = 1 > ' ,'=1J:=1 

i 2 l 2 

+ Mil "10 + "loll-ftl + ^2 II "20 + u 2o\\ n l 

Since 7iJ = 7iJ 0 © W}^, «i can be represented as u\ = 

uio + w^o where u 10 € Wj 0 and u^ 0 e H\q are orthogonal. Therefore, the minimizer to (28) must belong to the subspaces 
Therefore, H\ 0 zn&H\ Q . 

222 ^ pl u 88 m 8 (29) into (28) and setting the gradient with 

II wio + "io II = ll M ioll +ll"ioll ■ respect to c to 0, we see that c is the solution to (30). □ 
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